Method and tool for planning and dimensioning subsea pipeline-based transport systems for multiphase flows

ABSTRACT

This invention relates to a computer-implemented method for predicting fluid behaviour in pipeline-based multiphase flows, wherein the method comprises applying a one-dimensional computational fluid dynamic applying a finite volume method in the solver and which estimates the mass flux out of the finite control volumes by i) applying a polynomial to spatially reconstruct the mass present in each finite control volume, ii) reconstructing the flow velocity as a function of the x-component of the flow velocity vector to determine a domain of dependence for each finite control volume representing the distance the fluid has travelled during a time step, and iii) sum the spatially reconstructed mass being present in the domain of dependence for each finite control volume and assume the summarised mass passes out of the respective finite control volume over the applied time step.

CROSS REFERENCE TO RELATED APPLICATIONS

This application is a 371 National phase filing of PCT/EP2021/074668, filed Sep. 8, 2021, which application claims priority to Norwegian Patent Application No. 20201002, filed Sep. 11, 2020, the contents of which are hereby incorporated by reference in their entireties for all purposes.

FIELD OF INVENTION

This invention relates to a computer-implemented method for predicting fluid behaviour in pipeline-based transport systems for transport of multiphase flows.

BACKGROUND

New oil and gas reserves are often discovered a considerable distance apart from existing oil and gas extraction facilities leading to a need for transporting unprocessed fluids over long distances, maybe hundred kilometres or more. Multiphase fluid transport over such distances are often challenging due to complex non-linear interactions between the fluid phases resulting in an irregular and unstable behaviour causing pressure drops, deposit formations, and/or liquid accumulations which may manifest themselves as a range of problematic flow phenomena such as slug flow, bubbly flow, stratified flow, annular flow, chum flow, etc.

For example, liquid slugs are often problematic for long travelling distances such that being able to correctly predict the characteristics of the multiphase fluid flow, especially slug flow, is of great importance both in the design and operation of fluid transportation facilities of an oil field. E.g., the frequency and size of slugs are important parameters for the design of the fluid receiving facilities and can give important input to mechanical structural analysis securing the structural integrity. For gas-condensate type of fields this is involves equipment such as slug catcher and MEG-handling facilities. Here, underdesign may lead to severe operational problems shortening the lifetime of the facilities or production at a lower than planned plateau. Slugs are also problematic for the pipelines by causing load variations and subsequent vibrations that reduce the lifetime of pipe fittings and other vulnerable components. The ability to correctly predict the slug flow characteristics is therefore important both when designing the fluid transportation facilities and during operation of the fluid transportation facilities. In the latter case, the ability to predict the slug characteristics is important for investigating the effect of possible mitigating actions, such as e.g. topside choking, gas lift etc. Likewise, the ability to correctly predict liquid hold up is important. Underpredicting the pressure drop may lead to an undersized pipeline diameter and thus too small capacity during operation, while an overprediction of the pressure drop may lead to an oversized pipeline diameter which may create flow instabilities.

The ability to correctly predict multiphase fluid flow in transportation systems, especially when long travelling distances are involved, is vital for the economy and operability for e.g. hydrocarbon productions facilities. Accurate simulation models enabling predicting fluid flow behaviour consequences of different designs and approaches of the transportation system may save considerable investment and operation costs.

PRIOR ART

Prediction of the characteristics of multiphase flows requires specific computerbased simulations known as computational fluid dynamics (CFD). The computer codes of CFD-software usually consist of three main elements: (i) a pre-processor, (ii) a solver, and (iii) a post-processor [Ref. 1, pp. 1-4].

The pre-processor element concerns the definition/input of the fluid flow problem to be simulated by the operator via an operator interface and the transformation of the definition/input to a computer readable form applicable for the solver. The preprocessor element typically comprises:

-   -   input/definition of the geometry of the flow region of interest,         often denoted as the computational domain,     -   division of the computational domain into a set of         non-overlapping sub-domains, often denoted as grid generation,     -   input/defining the physical and chemical phenomena to be         simulated,     -   input/definition of the fluid properties, and     -   input/specification of boundary conditions at cells coinciding         with the boundary of the computational domain.

The post-processor element concerns the output of the simulation/simulation results and may include data-base modules for storing data, data visualisation modules, graphic presentation modules etc.

The solver element concerns the numerical solution of the natural laws/equations governing transport phenomena, convection, diffusion and source terms (creation of destruction of an entity of the flow. The numerical algorithm for solving the equations typically consists of three steps:

-   -   i) Integration of the governing equations of the fluid flow over         the computational domain.     -   ii) Discretisation—converting the integral equations into a         system of algebraic equations, and     -   iii) Solving the algebraic equations by iteration.

The governing equations of the fluid flow are mathematical statements expressing the conservation laws of physics to ensure conservation of mass, momentum and energy of the fluid. The fluid is treated as a continuum where its behaviour is described in terms of macroscopic properties such as velocity, pressure, density and temperature; see e.g. [Ref 1, pp. 9-10]. The conservation laws of physics lead to a transport equation describing the motion of an arbitrary scalar fluid property, p, which in a general formulation may be given on the form:

$\begin{matrix} {{\frac{\partial({\rho\phi})}{\partial t} + {{div}\left( {{\rho\phi}u} \right)}} = {{{div}\left( {D{{grad}\phi}} \right)} + S_{\phi}}} & (1) \end{matrix}$

where “ρ” is density, “u” is the velocity vector, “D” is the diffusion coefficient, “S_(ϕ)” is a source term for flow property ϕ, and “t” is the time coordinate. The first term (from left) describes the rate of increase of property ϕ of the fluid element, the second term describes the net rate of flow of fluid property ϕ out of the fluid element, the third term describes the net rate of increase of fluid property ϕ of the fluid element due to diffusion, and the fourth term describes the rate of increase of quantity ϕ due to sources.

The general formulated transport equation, eqn. (1), is usually a starting point for computational procedures in computational fluid dynamics by being applied to derive specific, often simplified, partial differential equations for the conservation of mass, momentum and energy for all fluid fields/phases of the fluid flow to be simulated in the specific geometry of the computational domain. CFD-simulations of gas-liquid flows in long pipelines usually requires using one-dimensional (1D) averaged conservation equations to achieve acceptable simulation times. For example, the mass conservation of a fluid phase, k, of a one-dimensional, unsteady and compressible multiphase flow with no source terms, eqn. (1) becomes reduced and may be given on the form:

$\begin{matrix} {{\frac{\partial{\overset{\hat{}}{\rho}(x)}_{k}}{\partial t} + \frac{{\partial{\overset{\hat{}}{\rho}(x)}_{k}}{u(x)}_{k}}{\partial x}} = 0} & (2) \end{matrix}$

Here “{circumflex over (ρ)}(x)_(k)” is the field mass of fluid phase k as a function of x, “u(x)_(k)” is the flow velocity for fluid phase k as a function of x, and {circumflex over (ρ)}_(k)=a_(k)ρ_(k) where a_(k) is the volume fraction of fluid phase k and ρ_(k) is the mass density of fluid phase k.

The computational domain specific derived conservation equations, such as e.g. eqn. (2), are integrated and discretised and solved over the computational domain by the solver element of the CFD-software. A widespread and much applied numerical solution technique is the finite volume method which is distinguished by treating the sub-domains of the computational domain as finite control volumes and integrating the governing equations of the fluid flow over each finite control volume (grid cell) of the computational domain. Applied on e.g. eqn. (2), the finite control volume integrates the eqn. over a finite control volume, ΔV_(i), and over a time step, Δt_(n), and may be given on the form

$\begin{matrix} {{{\int_{\Delta t_{n}}{{\frac{\partial}{\partial t}\left( {\int_{{\Delta V}_{i}}{{\overset{\hat{}}{\rho}(x)}_{k}dV}} \right)}dt}} + {\int_{\Delta t_{n}}{{\frac{\partial}{\partial t}\left( {\int_{{\Delta V}_{i}}{di{v\left( {{\overset{\hat{}}{\rho}(x)}_{k}{u(x)}_{k}} \right)}dV}} \right)}dt}}} = 0} & (3) \end{matrix}$

The next step in the finite volume method is transforming the continuous integral equation to a discrete algebraic equation. There are known several discretisation schemes for transforming the continuous integral equations into the discrete algebraic relations suitable for being solved numerically. The different difference schemes and their pros and cons in view of numerical stability and prediction accuracy are well known to the person skilled in the art. For flows dominated by convection, a well-known and much applied discretisation method for transforming the continuous integral equations into discrete algebraic relations is the upwind differencing scheme, see e.g. [Ref. 1, pp. 146-149], which applied on the mass conservation eqn. (2) may be given on the discrete form:

$\begin{matrix} {{\frac{\Delta{\overset{\hat{}}{\rho}}_{ki}}{\Delta t_{n}}\Delta V_{i}} = {F_{k,{i - {1/2}}} - F_{k,{i + {1/2}}}}} & (4) \end{matrix}$

Here, F_(k,i+1/2)={circumflex over (ρ)}_(k,i+1/2)u_(k,i+1/2)A_(i+1/2) is the mass flux on the right internal face of the i'th finite control volume, A_(i+1/2) is the cross-sectional area of the right internal face of the i'th finite control volume, ρ_(k,i+1/2) is the field mass of fluid phase k at position x_(i++1/2), and u_(k,i+1/2) is the fluid flow velocity at position x_(i+1/2). Depending on the choice of grid (e.g. collocated or staggered), the values at position x_(i+1/2) may have to be interpolated using methods well known to the person skilled in the art, for example using a first-order upwind scheme. Correspondingly, F_(k,i−1/2)={circumflex over (ρ)}_(k,i−1/2)u_(k,i−1/2)A_(i−1/2) is the mass flux on the left internal face of the i'th finite control volume. By “marching through” the entire computational domain, i.e. defining a discrete mass conservation equation for each of the entire set of i=1, . . . , N finite control volumes, the continuity equation for the conservation of mass is transformed to a set of algebraic relations of discrete values which may be solved by an explicit (direct) or implicit (indirect, i.e. iterative) solution algorithm with given boundary conditions and initial flow parameters.

Implicit numerical solution schemes determine the solution for the state of the system at a given time-step forward in time by solving an equation involving both the current state of the system (which is known) and the state of the system forward in time (which is unknown). This requires use of a matrix or iterative solution which gives somewhat extra computation but has the advantage that the numerical solution scheme is unconditionally stable and thus enables employing long timesteps which saves computation. However, at the cost of less accurate predictions.

Explicit numerical solution schemes determine the solution for the state of the system at a given time-step forward in time directly from the current state of the system (which is known). This has the advantage that the solution is significantly more accurate but is encumbered by being conditionally stable. A necessary condition for the stability of an explicit numerical solution scheme is that the numerical domain of dependence bounds the physical domain of dependence, i.e. the “physical” velocity u(x) should be less than the “marching velocity” Δx/Δt of the numerical method. This condition is known as the Courant-Friedrichs-Lewy (CFL) condition, which for a one-dimensional first order upwind scheme may be given as:

$\begin{matrix} {{\Delta t_{n}} \leq \frac{\Delta x_{i}}{u\left( {x_{i + {1/2}},t_{n + 1}} \right)}} & (5) \end{matrix}$

Here Δt_(n) is the n'th time step increment in the calculation, Δx_(i) is the length of the i'th finite control volume and u(x_(i+1/2), t_(n+1)) is the velocity of the fluid property at the upwind boundary of the i'th finite control volume at the end of the n'th time step. The CFL-condition should hold for every finite control volume of the computational domain and for every time step applied in the numerical calculations.

The limitation of the CFL-condition may be presented graphically such as e.g. given in FIGS. 2 a ) and 2 b). Both figures illustrate a section of a computational domain including five neighbouring finite control volumes having nodes at positions x_(i−3), x₁₋₂, . . . x_(i+1), respectively. The mass inside each finite control volume at the beginning of a n'th time step is indicated graphically by the vertical height of the node points at the centre of each finite control volume, i.e. the total amount of mass present in the finite control volume at the beginning of the time step equals the area of a rectangle being defined by the west and east internal faces, the abscissa and a horizontal line (parallel with the abscissa) crossing through the node point such as illustrated in FIG. 2 a ) by the shaded area in the i−2'th finite control volume. In the case of a first order upwind scheme, the mass passing through e.g. the internal face at position x_(i+1/2) during time step Δt_(n) is estimated as the area of the rectangle below the nodal point having a width determined as the product u(x_(i+1/2), t_(n+1))·Δt_(n), where u(x_(i+1/2), t_(n+1)) is the x-component of the velocity vector at the end of the n'th time step. If this product is equal to or less than the spatial increment Δx_(i)(i.e. the width of the i'th finite control volume), the CFL-condition is fulfilled. The amount of mass estimated passing out of the i'th finite control volume (over the internal face at position x_(i+1/2)) during the n'th time step is less than the total amount of mass being present in the finite control volume.

However, if the product u(x_(i+1/2), t_(n+1))·Δt_(n) is larger than the spatial increment Δx_(i), a non-physical situation arises because the numerical calculations estimate that more mass will flow out of the finite control volume than it contains as illustrated in FIG. 2 b ), which leads to estimates including “negative” masses. Thus, a violation of the CFL-condition may cause the estimation of negative masses in the finite control volumes which may lead to numerical instabilities in the solver.

The CFL-condition may impose a rather harsh numerical load in CFD-calculations of multiphase flows involving rapidly moving phases and/or flow phenomena requiring a fine-masked grid to be predicted acceptably accurate. Then the CFL-condition may dictate use of very small time-increments leading to unacceptably extensive and heavy computational loads to make CFD-simulations of multiphase fluid flows a practically applicable tool for designing and/or trouble-shooting transport systems for multi-phase flows.

Griebel and Klitz [Ref. 4] discloses a fast and mass-conserving method for simulation of two-phase flow problems. One popular method is the coupling of level-set and volume-of-fluid (CLSVOF), which benefits from the advantages of both approaches and results in improved mass conservation while retaining the straightforward computation of the curvature and the surface normal. Despite its popularity, details on the involved complex computational algorithms are hard to find and if found, they are mostly fragmented and inaccurate. In contrast, this article can be used as a comprehensive guide for an implementation of CLSVOF into the existing level-set Navier-Stokes solvers on Cartesian grids in three dimensions.

US 2019/228124 discloses an improved computer implemented method for modelling transport processes in fluids is disclosed. Instead of modelling based on using an infinitesimal fluid element of a continuous medium, the method approximates fluid flow in a fluid system as a model gas flow in a model gas system being identical to the fluid system. The method is adapted to model gas flow including dilute gas flow for high Knudsen numbers (Kn). The method delivers a new basis for prediction of dynamic evolution of the model gas system by considering a pre-established or known dynamic history of the system during a pre-initial period. A new generation of Computational Fluid Dynamics software products, which are based on the disclosed analytical tools and methods, are anticipated having capability to modelling gases from the continuum flow regime (Kn<0.01) to the free molecular flow regime (Kn>10), considerably higher accuracy of prediction, and lower computation cost.

From WO 2017/174532 it is known a simulation method (100) for simulating the thermo-fluid dynamic behavior of multiphase fluids in a hydrocarbons production and transport system, said method comprising the following steps: outlining (110) the hydrocarbons production and transport system as a plurality of interconnected component blocks, thus creating a schematic representation; modelling (120) each component block with a simplified analytical mathematical model selected from the group of models comprising at least one conduit model, a valve model, a reservoir model and a separator model, each simplified analytical mathematical model comprising a plurality of constitutive equations adapted to describe the thermo-fluid dynamic behaviour of the corresponding component block; generating (130) an oriented graph on the basis of the schematic representation; determining (140) a plurality of topological equations on the basis of the oriented graph; determining (150) a plurality of output variables adapted to describe the thermo-fluid dynamic behaviour of the system by solving the set of the plurality of topological equations and of the constitutive equations.

WO 2019/164859 discloses a compositional simulator that is fully compositional and more robust and accurate is disclosed. Relative permeability (kr) and capillary pressure (Pc) are modelled as state functions, making them unique for a given set of inputs, which can include Euler characteristic, wettability, pore connectivity, saturation, and capillary number. All of these are made to be a function of composition, T, and P or rock properties. These state function kr-Pc models are fully compositional and can fit experimental data, including complex processes such as hysteresis. The models can be tuned to measured relative permeability data, and then give consistent predictions away from that measured data set. Phase labelling problems are eliminated. Flux calculations from one grid block to another are based on phase compositions. Simulations for three or four-phase hydrocarbon phases are possible. Time-step sizes increase to stability limits of implicit-explicit methods. Fully implicit methods are possible and significant improvements are expected.

From US 2018/260499 it is known a fluid is modelled as a set of discrete particles. Each of the particles is associated with one or more properties, and a spatial distance comprising a smoothing length over which the one or more properties are to be smoothed. A corresponding trajectory is simulated for each of the particles. The corresponding trajectory is used to formulate a first solution for simulating transport within the fluid. A first predicted error is determined for the first solution. An iterative adjustment is performed to at least one of: a quantity of particles, the smoothing length, or the one or more corresponding properties, to formulate a second solution for simulating transport with the fluid, and a second predicted error is determined for the second solution, until the second predicted error is within a predefined boundary.

OBJECTIVE OF THE INVENTION

The main objective of the invention is to provide a computer implemented method for predicting fluid behaviour of a multiphase flow in a pipeline-based transport system having a solver enabling using an explicit numerical solution scheme which is numerically stable independently of the Courant-Friedrichs-Lewy (CFL) condition when predicting the mass flow.

A further objective of the invention is to provide a computer implemented method for predicting fluid behaviour of a multiphase flow in a pipeline-based transport system which preserves the positivity of the mass everywhere.

Another objective of the invention is to provide a computer implemented method for designing transport systems for multiphase fluid flows.

A further objective of the invention is a computer implemented simulation tool for designing/optimising and/or trouble-shooting a pipeline-based transport system for multi-phase fluid flows.

DESCRIPTION OF THE INVENTION

The present invention is based on the realisation that the positivity of the mass for one dimensional flows may be preserved in the numerical calculations independently of the Courant-Friedrichs-Lewy (CFL) condition by applying an approach for estimating the mass flux which comprises for a n'th time step, Δt_(n);

-   -   applying a polynomial to spatially reconstruct the mass,         {circumflex over (ρ)}_(k)(x), for the k'th fluid phase being         present in each of the N finite control volumes of the         computational domain, and then     -   applying the discretised flow velocity, u_(k) (which can be         located either at the centre point of the mass cells in a         collocated grid, or at the internal faces in a staggered grid),         to determine a domain of dependence,         _(k,j), representing the distance the k'th fluid phase has         travelled from its starting point to a j'th internal face, where         jϵ1/2, 3/2, . . . , N−1/2, during the n'th time step for each         internal face j of the computational domain, and     -   sum the spatially reconstructed mass being present in the domain         of dependence,         _(k,j) and assume it passes through the j'th internal face         during the n'th time step.

As used herein, the term “computational domain” is a virtual representation of the geometry of a section of a pipe containing a multiphase flow which is to be simulated. The computational domain is divided into a number of N non-overlapping finite control volumes where each i'th, i=1, . . . , N, finite control volume represents a separate “slice” of the computational domain. As used herein, index “i” is applied when the focus is on the finite control volumes (iϵ1, . . . , N) and index “j” is applied when the focus is on the internal faces jϵ1/2, . . . , N+1/2. The notations “i−1/2” and “i+1/2” as used herein relate to the internal faces of the finite control volume i (left and right internal faces, respectively), and “j−1/2” and “j+1/2” are used to relate finite control volumes to the internal face j (left and right control volumes, respectively). The terms “left” and “west” are used herein interchangeably to indicate a direction opposite to the grid marching direction (i.e. towards decreasing index “i” or “j”), and consequently, the terms “right” and “east” are used herein interchangeably to indicate a direction in the grid marching direction (i.e. towards increasing index “i” or “j”).

The spatial reconstruction of the mass of the k'th fluid phase over the computational domain may be obtained by any manner known and conceivable to the skilled person, provided the following conditions are satisfied:

-   -   The reconstructed mass function should have, within each finite         control volume, an average value equal to the discretised mass         at the cell's centre point. In other words, it should conserve         the mass present in each cell

In one example embodiment, the spatial reconstruction of the mass of the k'th fluid phase over the computational domain may also satisfy the condition that there must not be any negative value anywhere in the reconstruction.

FIG. 3 a ) illustrates schematically an example of using a first order polynomial to represent the mass present in the finite control volumes as linear interpolation lines over each finite volume. Use of higher order polynomials gives a more accurate spatial reconstruction.

The range and length of the domain of dependence,

_(k,j), depends on the flow direction, the flow velocity, and the duration of the n'th time step. If the flow direction is positive, the domain of dependence stretches from the j'th internal face in a westward direction as shown in FIGS. 3 b ) and 3 c), first in cell i (the j'th internal face is between control volumes i and i+1). As seen on the figures, the domain of dependence,

_(k,j), may be longer than the i'th finite control volume such that it extends into one or more neighbouring finite control volumes. In the example shown schematically in FIGS. 3 a ) and 3 b), the domain of dependence stretches into two and a half neighbouring finite control volumes from position x_(0,k) ^(n) being inside the i−3'th finite control volume to position x_(j). If the flow direction is negative, the domain of dependence stretches in eastward direction (not shown on the figures).

FIG. 3 c ) illustrates the mass being present in the domain of dependence,

_(k,j), by the shaded area between the abscissa and the graphic representation of {circumflex over (ρ)}_(k)(x) (in this example, linear interpolation lines). The shaded area over the domain of dependence represents the sum of mass which is assumed passing through the j'th internal face during the n'th time step. Since all mass present in the domain of dependence is summed and assumed passing through the corresponding internal face, the above approach ensures that the mass passing through the internal face actually exists, thus not pulling more mass out of a control volume than is present. Furthermore, since the domain of dependence may extend over more than one finite control volume, the numerical method is stable also when the “physical” velocity, u(x) becomes larger than the “marching velocity”, Δx/Δt, of the numerical scheme.

Thus, in a first aspect, the invention relates to a computer implemented method for predicting fluid behaviour of a multiphase flow in a pipeline-based transport system where the flow contains a number of k, where k is a positive integer, fluid phases, wherein the method comprises:

-   -   applying a one-dimensional (1D) computational fluid dynamic         (CFD) model describing the geometry of a section of interest of         the pipeline-based transport system and the multiphase flow         flowing therein,     -   solving the 1D CFD model to simulate the fluid behaviour of the         multiphase flow in the section of interest of the pipeline-based         transport system, and describing the determined fluid behaviour         by presenting the simulation results as macroscopic fluid         properties such as flow velocity, pressure, density,         temperature, etc., and     -   the 1D CFD model applies a finite volume method to solve the         model, wherein the geometry of the section of interest of the         pipeline-based transport system is defined as a computational         domain extending along an axis represented by the cartesian         coordinate x and being divided into a set of N, where N is a         positive integer, non-overlapping finite control volumes         separated by an internal face between adjacent finite control         volumes, characterised in that     -   the one-dimensional computational fluid dynamic model is adapted         to estimate the mass flow of a k'th fluid phase out of a i'th         finite control volume during a n'th time step by applying a         polynomial to spatially reconstruct the mass, {circumflex over         (ρ)}_(k)(x), of the k'th fluid phase being present in each of         the N finite control volumes of the computational domain, and         then     -   for each j'th internal face, where jϵ1/2, . . . , N+1/2, of the         computational domain, execute the following steps i) and ii):     -   i) reconstruct the flow velocity, u_(k)(x), of the k'th fluid         phase as a function of position x and apply the reconstruction         to determine a domain of dependence,         _(k,j), representing the distance the k'th fluid phase has         travelled during the n'th time step, and     -   ii) sum the spatially reconstructed mass being present in the         domain of dependence,         _(k,j), and assume the summarised mass passes through the j'th         internal face over the n'th time step, into the i'th finite         control volume when u_(k)(x_(j))<0 or into the i+1'th finite         control volume when u_(k)(x_(j))>0, where u_(k)(x_(j)) is the         flow velocity at the j'th internal face.

As used herein, the subscript n denoting the n'th time step is omitted to simplify the notation since all relations relate to the n'th time step. The concepts of the finite volume method, explicit or implicit numerical solution algorithms, upwind differencing scheme, staggered or collocated grids, etc. and how to implement these in CFD-models are well known and mastered by persons skilled in the art.

The present invention is not limited to any choice of discretisation scheme or numerical solution algorithm except that the CFD-model shall apply a finite volume method approach. That is, the numerical scheme according to the first aspect of the invention may be applied in any 1D CFD-model regardless of which differencing scheme and/or explicit or implicit numerical solution algorithm being implemented in the CFD-model. However, the numerical scheme according to the invention is especially beneficial for explicit numerical schemes by enabling violating the CFL-condition without excessively compromising the numerical stability, making CFD-models having an explicit numerical scheme a preferred example embodiment. Explicit formulated numerical schemes have a relatively high numerical accuracy and are well suited for multicore CPUs and computational cluster machines frequently used in cloud computing.

The numerical scheme according to the first aspect of the invention relates to the solution of the transport equations for the mass transfer of each fluid phase and field present in the multiphase flow. The transport equations governing momentum and energy of the fluid phases and fields of the multiphase flow may be solved by the CFD-model in any known manner. Thus, the present invention encompasses any CFD-model having a solver for solving the transport equations where the transport equations governing the mass conservation/flow of mass is solved according to the first aspect of the invention, and where the transport equations governing the momentum and energy are solved in any known manner. The transport equations governing momentum and energy may be subject to a CFL-condition.

As used herein, the term “pipeline-based transport system” encompasses all components of the transport system necessary to transport the fluid including pipeline segments, splits, joins, valves, pumps, sources, sinks, etc. An example of a pipeline-based transport system for produced liquids in oil and gas extraction is shown in FIG. 1 which illustrates schematically an example embodiment of such transportation system. This example embodiment comprises a plurality of tiebacks/pipelines (2) connecting a production well (1) to a nearby satellite hub (3) which collects the produced fluid in a region and passes the produced fluid in a satellite pipeline (4) to a common hub (5). The example embodiment comprises four satellite hubs (3) connected to the common (5) by a satellite pipeline (4) each. The common hub (5) passes the produced fluid to a processing facility located either offshore on the sea surface via a riser (not shown in this embodiment) or to an onshore production facility via fluid transportation pipeline (6). The transport system usually involves one or more fluid pumps (7) to provide the necessary flow pressure to move the fluids through the transport system.

The reconstruction of the flow velocity, u_(k)(x), of the k'th fluid phase as a function of position x to determine a domain of dependence,

_(k,j), may be obtained in several ways known to the skilled person. The invention according to the first aspect may apply any conceivable manner known to the skilled person.

As an example, linear interpolation of the discretised flow velocity within each finite control volume may be applied to reconstruct the flow velocity, u_(k)(x), as a function of the position x, here given for the i'th finite control volume and a staggered mesh arrangement:

$\begin{matrix} {{{u(x)} = {{\frac{x_{i + {1/2}} - x}{x_{i + {1/2}} - x_{i - {1/2}}}u_{k,{i - {1/2}}}} + {\frac{x - x_{i - {1/2}}}{x_{i + {1/2}} - x_{i - {1/2}}}u_{k,{i + {1/2}}}}}},} & (6) \end{matrix}$ x ∈ [x_(i − 1/2), x_(i + 1/2)]

Eqn. (6) may be rearranged to the form:

$\begin{matrix} {{{u_{k}(x)} = {{\overset{\_}{u}}_{k.i} + {\frac{{\Delta u}_{k,i}}{{\Delta x}_{i}}\left( {x - {\overset{\_}{x}}_{i}} \right)}}},} & (7) \end{matrix}$ x ∈ [x_(i − 1/2), x_(i + 1/2)] Here ${{\overset{\_}{x}}_{i} = \frac{x_{i + {1/2}} + x_{i - {1/2}}}{2}},$ ${{\overset{\_}{u}}_{k,i} = \frac{u_{k,{i - {1/2}}} + u_{k,{i + {1/2}}}}{2}},$ Δx_(i) = x_(i + 1/2) − x_(i − 1/2), and Δu_(k, i) = u_(k, i + 1/2) − u_(k, i − 1/2).

Equations (6) and (7) can be written for all the finite control volumes. For example, if they are written for the i+1'th finite control volume, they will be valid in x∈[x_(i+1/2), x_(i+3/2)]. By joining all the finite control volumes, a reconstruction for any x in the computational domain is obtained.

The domain of dependence,

_(k,j), for the j'th internal face at the n'th time step may be determined by applying and integrating the equation of the flow characteristics:

$\begin{matrix} {{\frac{dx}{dt} = {u_{k}(x)}},} & (8) \end{matrix}$

with the velocity as expressed in eqn. (7). This gives the following relation:

$\begin{matrix} {{x_{k}(t)} = {{\overset{\_}{x}}_{i} + {e^{- \frac{{\Delta u}_{k,i}}{{\Delta x}_{i}}}\left( {x_{s,k}^{n} - {\overset{\_}{x}}_{i}} \right)} + {{\overset{\_}{u}}_{k,i}\frac{{\Delta x}_{i}}{{\Delta u}_{i}}\left( {e^{- \frac{{\Delta u}_{k,i}}{{\Delta x}_{i}}} - 1} \right)}}} & (9) \end{matrix}$

where x_(s,k) ^(n) is the starting position of the flow characteristic of the k'th fluid phase at the beginning of the applied time step, crossing the position x_(k)(t) at a time t. The finite control volume index i denotes the cell in the upwind direction and is equal to j−1/2 when the flow direction is positive, i.e. when u_(k)(x_(j))>0. When the flow direction is negative, i.e. when u_(k)(x_(j))<0, eqn. (11) still applies if the index i is set to be j+1/2. Eqn. (9) may be solved for x_(s,k) ^(n), under the requirement that x_(k)(Δt)=x_(j), which means that at the end of the time step Δt, the characteristic will cross the internal face j. This gives:

$\begin{matrix} {x_{s,k}^{n} = {{\overset{\_}{x}}_{i} + {e^{{- \frac{\Delta u_{k,i}}{\Delta x_{i}}}\Delta t}\left( {x_{j} - {\overset{\_}{x}}_{i}} \right)} + {{\overset{\_}{u}}_{k,i}\frac{\Delta x_{i}}{\Delta u_{i}}\left( {e^{{- \frac{\Delta u_{k,i}}{\Delta x_{i}}}\Delta t} - 1} \right)}}} & (10) \end{matrix}$

If the starting position x_(s,k) ^(n) lies within the i'th finite volume, the domain of dependence,

_(k,j), for the k'th fluid phase and the j'th internal face at the n'th time step, will extend from the starting position x_(s,k) ^(n) to the j'th internal face. In this case it may be advantageous to denote x_(s,k) ^(n) as x_(0,k) ^(n) to indicate that it is the starting point of the domain of dependence,

_(k,j):

k , j = { [ x 0 , k n , x j , ] if ⁢ u k , j > 0 [ x j , x 0 , k n ] if ⁢ ⁢ u k , j < 0 ( 11 )

One advantage of the invention according to the first aspect is that it enables applying time steps being larger than the time the flow characteristic uses to cross the i'th finite control volume without compromising the numerical stability. However, if the n'th time step is larger than the time the flow characteristic uses to cross the i'th finite control volume, the starting position determined by eqn. (10) will be lying outside of the i'th finite control volume. To accommodate for this situation, the invention according to the first aspect may further comprise a procedure to determine the domain of dependence which takes into account the time needed for the k'th flow characteristic to cross the finite control volumes. The time needed for the flow characteristic to cross a finite control volume may be determined by solving eqn. (10) for Δt, which when applied in the i'th finite control volume may be given as:

$\begin{matrix} {{\Delta t_{c,k,i}} = {\frac{\Delta{x}_{i}}{\Delta u_{k,i}}{\ln\left( \frac{u_{k,{i + {1/2}}}}{u_{k,{i - {1/2}}}} \right)}}} & (12) \end{matrix}$

Here Δt_(c,k,i) is the time needed for the k'th flow characteristic to cross the i'th finite control volume. Eqn. (12) may also be applied to determine the time needed for the k'th flow characteristic to cross the i−1'th finite control volume by applying the variables related to the i−1'th finite control volume: Δx_(i−1), Δu_(k,i−1), u_(k,i−1/2) and u_(k,i−3/2). Eqn. (12) may be applied similarly to determine the crossing time for any other finite control volume. If the ratio of velocities is negative, Eqn. (12) becomes invalid. It happens when the velocity changes sign in the cell. In this case, no characteristic can cross the position where u_(k)(x)=0. Thus, the crossing time is infinite, and the characteristic will start within the cell.

Close to the boundary nodes the domain of dependence may extend out of the computational domain

. This means that the total crossing time from the computational domains' west or east end to the internal face j is smaller than the time step Δt. Then, there is a fraction of the time step for which the mass crossing internal face j will originate from outside the computational domain. To accomplish for this situation, the method may be extended with a special treatment of the nodes. Depending on the type of boundary condition that we want to achieve, two quantities may be useful to calculate, the remainder of the time step Δt_(r,k), and the extension of the domain of dependence outside of the computational domain.

When the domain of dependence

_(k,j) would stretch outside of the computational domain,

, the remainder of the time step Δt_(r,k) is defined as

$\begin{matrix} {{\Delta t_{r,k}} = {{\Delta t} - \left\{ \begin{matrix} {{\sum}_{i \in {\lbrack{{j + 1},N}\rbrack}}\Delta t_{c,k,i}} & {{{if}u_{k,j}} < 0} \\ {{\sum}_{i \in {\lbrack{1,j}\rbrack}}\Delta t_{c,k,i}} & {{{if}u_{k,j}} > 0} \end{matrix} \right.}} & \left( {12a} \right) \end{matrix}$

For the extension of the domain outside of the computational domain,

, the velocity outside of the pipe must be extrapolated. As an example, it can be constantly extrapolated and equal to the velocity in the west- or eastmost cell, for the west or east node, respectively. Then, equation (10) can be applied with the remainder of the time step Δt_(r,k), and ū_(k,j) set equal to the extrapolated velocity. Equation (10) in this case will have a numerical issue, since there is Δu_(k,i)=0 at the denominator due to the constant extrapolation. This issue is solved further down.

The notation for intervals as used herein follows the international standard ISO 80000-2, where the brackets “[” and “]” indicate a closed interval border, and the parenthesises “(” and “)” indicate an open interval border. For example, [a, b] is the closed interval containing every real number from a included to b included: [a, b]={x∈

|a≤x≤b}, while (a, b] is the left half-open interval from a excluded to b included: (a,b]={x∈

|a<x≤b}.

An example embodiment of a procedure for determining the domain of dependence of the j'th internal face and k'th fluid phase when u_(k,j)≠0, may be given as:

Initialisation:

-   -   i) if u_(k,j)>0:         -   set the upwind cell index i=j−1/2 and the direction index             s=−1     -   if u_(k,j)<0:         -   set the upwind cell index i=j+1/2 and the direction index             s=1     -   (ii) set the cell counter m=0     -   (iii) set Δt_(r), the remainder of the time step after crossing         cells i to i+sm, equal to the full time step Δt

Step 1:

-   -   (iv) apply eqn. (12) to determine the crossing time of the         (i+sm)'th finite control volume, Δt_(c,k,i+sm)     -   (v) if Δt_(c,k,i+sm)>Δt_(r) (the characteristic will not cross         the whole (i+sm)'th finite control volume:         -   Go to step 2     -   or else (the characteristic will cross the whole (i+sm)'th         finite control volume and enter next one):         -   if i+s(m+1)<1 or i+s(m+1)>N: (the end of the pipe's domain             has been reached)             -   set Δt_(r)=Δt_(r)−Δt_(c,k,i+sm), and terminate the                 procedure,         -   or else:             -   set m=m+1 and Δt_(r)=Δt_(r)−Δt_(c,k,i+sm), and go to                 step 1

Step 2:

-   -   (vi) Apply eqn. (10) with the flow velocity shorthand's         ū_(k,i+sm) and Δu_(k,i+sm) as defined below eqn. (7), and         timestep Δt=Δt_(r) to determine the characteristic starting         position, x_(0,k) ^(n), which defines the domain of dependence         _(k,j) according to eqn. (11) and terminate the procedure.

The above procedure for determining the domain of dependence,

_(k,j), may be presented graphically as shown in e.g. FIGS. 3 b ) and 3 c), which illustrates a domain of dependence stretching from the j'th internal face across three neighbouring finite control volumes and a distance into a fourth neighbouring finite control volume. The flow direction is indicated by the black arrows located at the internal faces. The figure illustrates an example with a positive flow direction such that the domain of dependence stretches from the internal face at position x_(j) to the starting position x_(0,k) ^(n) lying inside the i−3'th finite control volume.

The above procedure for determining the domain of dependence is described for either a negative or a positive flow direction. In the case of zero flow velocity there is no flow to be modelled and the domain of dependence will be zero. This situation may be taken care of by setting the domain of dependence

_(k,j), to zero when u_(k,j)=0. The particular case where the velocity change sign from one internal face to the next is naturally handled, since the condition to leave step 1, Δt_(c,k,i+sm)>Δt_(r), where Δt_(c,k,i+sm)=∞, will be true. Step 2 has no particular issue with velocity changing sign.

The above procedure for determination of the domain of dependence may have an issue when the velocities in cells i−1/2 and i+1/2 are equal to each other, since the denominator Δu_(k,j) of eqns. (10) and (12) goes to towards zero. However, taking the limit of constant velocity we have that:

$\begin{matrix} {{\lim\limits_{{\Delta u}\rightarrow 0}{x(t)}} = {x_{0,k}^{n} + {{\overset{\_}{u}}_{k,i}t}}} & (13) \end{matrix}$

which shows that the Δu_(k,i) in the denominator is a numerical issue rather than a fundamental problem with the above approach. The above approach may thus be made numerically robust by e.g. defining the change in the Courant number over the i'th finite control volume, ΔC_(i), to be:

$\begin{matrix} {{{\Delta C_{k,i}} = {\Delta u_{k,i}\frac{\Delta t}{\Delta x_{i}}}},} & (14) \end{matrix}$

and apply the expression of eqn. (14) to rewrite the problematic term in eqn. (10) as:

$\begin{matrix} {{\overset{\_}{u}}_{k,i}\frac{\left( {\left( e^{{- \Delta}C_{k,i}} \right) - 1} \right)}{\Delta C_{k,i}}\Delta t} & (15) \end{matrix}$

and expand the expression in a Taylor series:

$\begin{matrix} {\frac{\left( {1 - e^{{- \Delta}C_{k,i}}} \right)}{\Delta C_{k,i}} = {1 - {\frac{1}{2}\Delta C_{k,i}} + {\frac{1}{6}\Delta C_{k,i}^{2}} + {O\left( {\Delta C_{k,i}^{3}} \right)}}} & (16) \end{matrix}$

If the second order term in eqn. (16) is less than the ϵ value returned from the Fortran function EPSILON which returns a positive real value being the smallest possible value of ϵ making 1+ϵ not equal to 1 on the computer system, then the numerical value of the function on the left side of eqn. (16) is numerically indistinguishable from the two first terms of the Taylor expansion, such that a numerically robust form of eqn. (10) may be given as:

$\begin{matrix} {x_{s,k}^{n} = {{\overset{\_}{x}}_{i} + {e^{{- \Delta}C_{k,i}}\left( {x_{j} - {\overset{\_}{x}}_{i}} \right)} - {{\overset{\_}{u}}_{k,i}\Delta t\left\{ \begin{matrix} {1 - {\frac{1}{2}\Delta C_{k,i}}} & {{{if}{❘{\Delta C_{k.i}}❘}} < \sqrt{6\epsilon}} \\ \frac{\left( {1 - e^{{- \Delta}C_{k,i}}} \right)}{\Delta C_{k,i}} & {otherwise} \end{matrix} \right.}}} & (17) \end{matrix}$

where the index i=j−1/2 when the flow direction is positive, i.e. when u_(k)(x_(j))>0. When the flow direction is negative, i.e. when u_(k)(x_(j))<0, eqn. (17) still applies if the index i is set to be: i=j+1/2. At is the applied time step.

Similarly, the time needed for the flow characteristic to cross a finite control volume may be made numerically robust by defining r as the ratio:

$\begin{matrix} {{r_{k,i} = \frac{u_{k,{i + {1/2}}}}{u_{k,{i - {1/2}}}}},} & (18) \end{matrix}$

and rewrite eqn. (12) on the form:

$\begin{matrix} {{{\overset{\_}{u}}_{k,i}\frac{\Delta t_{c,1}}{\Delta x_{i}}} = {\frac{r_{k,i} + 1}{2}\frac{\ln\left( r_{k,i} \right)}{\left( {r_{k,i} - 1} \right)}}} & (19) \end{matrix}$

The Taylor expansion of the right-hand side of eqn. (19) around r_(k,i)=1 is:

$\begin{matrix} {{\frac{r_{k,i} + 1}{2}\frac{\ln\left( r_{k,i} \right)}{\left( {r_{k,i} - 1} \right)}} = {1 + {\frac{1}{12}\left( {r_{k,i} - 1} \right)^{2}} - {\frac{1}{12}\left( {r_{k,i} - 1} \right)^{3}} + {\frac{3}{40}\left( {r_{k,i} - 1} \right)^{4}} - {O\left( {r_{k,i} - 1} \right)}^{5}}} & (20) \end{matrix}$

If the second order term of eqn. (20) is less than the ϵ value returned from the Fortran function EPSILON, then the numerical value of the function on the left of eqn. (20) is indistinguishable from 1, such that a numerically robust expression for the flow characteristic to cross a finite control volume, Δt_(c,1), may be given as:

$\begin{matrix} {{{\overset{\_}{u}}_{k,i}\frac{\Delta t_{c,1}}{\Delta x_{i}}} = \left\{ \begin{matrix} 1 & {{{if}{❘{r_{k,i} - 1}❘}} < {2\sqrt{3\epsilon}}} \\ {\frac{r_{k,i} + 1}{2}\frac{\ln\left( r_{k,i} \right)}{\left( {r_{k,i} - 1} \right)}} & {otherwise} \end{matrix} \right.} & (21) \end{matrix}$

In an example embodiment, the invention according to the first aspect may further comprise determining the domain of dependence,

_(k,j) by executing the following steps i) to vi):

Step 1 (initialisation):

-   -   i) if u_(k,j)>0:         -   set an upwind cell index i=j−1/2 and a direction index s=−1,     -   if u_(k,j)<0:     -   set an upwind cell index i=j+1/2 and a direction index s=1,     -   ii) set a cell counter m=0,     -   iii) set Δt_(r), the remainder of the time step after crossing         cell i to i+sm, equal to the full time step Δt,

Step 2:

-   -   iv) if |r_(k,i+sm)−1|<2√{square root over (3ϵ)}:

${{{set}\Delta t_{c,k,{i + {sm}}}} = \frac{\Delta x_{i + {sm}}}{{\overset{\_}{u}}_{k,{i + {sm}}}}},$

-   -   -   otherwise:

${{{set}\Delta t_{c,k,{i + {sm}}}} = {\frac{\Delta x_{i + {sm}}}{{\overset{\_}{u}}_{k,{i + {sm}}}}\left( {\frac{r_{k,{i + {sm}}} + 1}{2}\frac{\ln\left( r_{k,{i + {sm}}} \right)}{\left( {r_{k,{i + {sm}}} - 1} \right)}} \right)}},$

-   -   -   where

${r_{k,{i + {sm}}} = \frac{u_{k,{i + {sm} + {1/2}}}}{u_{k,{i + {sm} - {1/2}}}}},$

-   -   -    and             -   ϵ is a positive real number obtained from executing a                 Fortran EPSILON computer implemented function,

    -   v) if Δt_(c,k,i+sm)>Δt_(r):         -   go to step vi),

    -   or else if i+s(m+1)<1 or i+s(m+1)>N:         -   set Δt_(r)=Δt_(r)−Δt_(c,k,i+sm), and terminate the             procedure,

    -   or else:         -   set m=m+1 and Δt_(r)=Δt_(r)−Δt_(c,k,i+sm), and go to step             iv),

Step 3):

-   -   vi) if |C_(k,i+sm)|<√{square root over (6ϵ)}:         -   determine a characteristic starting position, x_(0,k) ^(n),             by solving the following eqn;

${x_{0,k}^{n} = {{\overset{\_}{x}}_{i + {sm}} + {e^{{- \Delta}C_{k,{i + {sm}}}}\left( {x_{j + {sm}} - {\overset{\_}{x}}_{i + {sm}}} \right)} - {{\overset{\_}{u}}_{k,{i + {sm}}}\Delta{t_{r}\left( {1 - {\frac{1}{2}\Delta C_{k,{i + {sm}}}}} \right)}}}},$

-   -   -    and         -   go to step vii),

    -   Otherwise:         -   determine a characteristic starting position, x_(0,k) ^(n),             by solving the following eqn.;

${x_{0,k}^{n} = {{\overset{¯}{x}}_{i + {sm}} + {e^{{- \Delta}C_{k,{i + {sm}}}}\left( {x_{j + {sm}} - {\overset{¯}{x}}_{i + {sm}}} \right)} - {{\overset{¯}{u}}_{k,{i + {sm}}}\Delta{t_{r}\left( \frac{\left( {1 - e^{{- \Delta}C_{k,{i + {sm}}}}} \right)}{\Delta C_{k,{i + {sm}}}} \right)}}}},$

-   -   -    and         -   go to step vii),

    -   where

${{\overset{¯}{x}}_{i + {sm}} = \frac{x_{i + {sm} + {1/2}} + x_{i + {sm} - {1/2}}}{2}},$ ${{\overset{¯}{u}}_{k,{i + {sm}}} = \frac{u_{k,{i + {sm} - {1/2}}} + u_{k,{i + {sm} + {1/2}}}}{2}},$ Δx_(i + sm) = x_(i + sm + 1/2) − x_(i + sm − 1/2), Δu_(k, i + sm) = u_(k, i + sm + 1/2) − u_(k, i + sm − 1/2), ${{\Delta C_{k,{i + {sm}}}} = {\Delta u_{k,{i + {sm}}}\frac{\Delta t}{\Delta x_{i + {sm}}}}},$

-   -   -    and

Step 4:

-   -   vii) if u_(k,j)>0:         -   set             _(k,j)=[x_(0,k) ^(n),x_(j)],     -   or if u_(k,i)<0:         -   set             _(k,j)=[x_(j), x_(0,k) ^(n)].

This terminates the procedure.

As mentioned above, the spatial reconstruction of the mass in the finite control volumes may be obtained in several ways known to the person skilled in the art. The invention according to the first aspect may apply any mathematical technique for spatially reconstructing the mass in the finite control volumes known and conceivable to the skilled person, provided the following condition is satisfied:

-   -   The reconstructed mass function should have, within each finite         control volume, an average value equal to the discretised mass         at the cell's centre point. In other words, it should conserve         the mass present in each cell

In one example embodiment, the spatial reconstruction of the mass of the k'th fluid phase over the computational domain may also satisfy the condition that there must not be any negative value anywhere in the reconstruction.

The probably simplest spatial reconstruction is a piecewise constant function, equal within each finite control volume to the discretised mass at the finite control volume centre point. Another example is using a piecewise linear function composed of first order polynomial, one per finite control volume, whose average over the volume is equal to the discretised mass at the finite control volume centre point. Higher order polynomials can also be used and will result in a higher convergence order of the scheme. Especially even-order polynomials are suited since they allow expanding the stencil equally to both sides of the finite control volume being subject for mass reconstruction. Thus, the spatial reconstruction of the mass as a function of position variable x may be obtained by applying an even numbered polynomial:

{circumflex over (ρ)}_(k,i)*(x)=Σ_(α=0) ^(β) c _(α) x ^(α),β∈[0,2,4, . . . ]  (22)

and then, for each i'th finite control volume, where iϵ1, . . . , N, of the computational domain, and each k'th phase, define a set of coefficients, c_(α), through the condition:

$\begin{matrix} {{\int_{x_{i - {1/2}}}^{x_{i + {1/2}}}{{{\overset{\hat{}}{\rho}}_{k,i}^{*}(x)}{dx}}} = {{\overset{\hat{}}{\rho}}_{k,i}\Delta x_{i}}} & (23) \end{matrix}$

and solve for the coefficients c₀, c₁, . . . , c_(β) to express the spatial reconstruction of the mass of k'th phase present in the i'th finite control volume as {circumflex over (p)}_(k,i)*(x)=c₀+c₁x+c₂x²+ . . . +c_(β)x^(β). There is one set of coefficients c₀, c₁, . . . , c_(β) for each phase k in each cell i, but the k and i indices are dropped for readability. By repeating this procedure for every finite control volume of the computational domain, the mass of each finite control volume is spatially reconstructed over the entire computational domain.

The simplest spatial reconstruction is to set β=0 such that p(x) becomes the average mass of phase k in the i'th finite control volume at the beginning of the n'th time step, {circumflex over (ρ)}_(k,i) ⁰.

An example embodiment employing a second order polynomial; p(x)=c₀+c₁x+c₂x² is as follows. The coefficients c₀, c₁, c₂ for phase k in the i'th cell are calculated by requiring that the integrals of the polynomial in the three ranges xϵ[x_(i−3/2), x_(i−1/2)], xϵ[x_(i−1/2), x_(i+1/2)] and xϵ[x_(i+1/2), x_(i+3/2)], be equal to the masses in the i−1'th, the i'th, and the i+1'th finite control volumes, respectively. Using that:

$\begin{matrix} {{{\int_{x_{i - {1/2}}}^{x_{i + {1/2}}}c_{0}} + {c_{1}x} + {c_{2}x^{2}{dx}}} = {{c_{0}\Delta x_{i}} + {c_{1}{\overset{¯}{x}}_{i}\Delta x_{i}} + {c_{2}\left( {{{\overset{\_}{x}}_{i}^{2}\Delta x_{i}} + {\frac{1}{12}\Delta x_{i}^{3}}} \right)}}} & (24) \end{matrix}$

the system to solve becomes:

$\begin{matrix} {{{c_{0}\Delta x_{i - 1}} + {c_{1}{\overset{¯}{x}}_{i - 1}\Delta x_{i - 1}} + {c_{2}\left( {{{\overset{¯}{x}}_{i - 1}^{2}\Delta x_{i - 1}} + {\frac{1}{12}\Delta x_{i - 1}^{3}}} \right)}} = {{\overset{\hat{}}{\rho}}_{k,{i - 1}}\Delta x_{i - 1}}} & (25) \end{matrix}$ $\begin{matrix} {{{c_{0}\Delta x_{i}} + {c_{1}{\overset{¯}{x}}_{i}\Delta x_{i}} + {c_{2}\left( {{{\overset{\_}{x}}_{i}^{2}\Delta x_{i}} + {\frac{1}{12}\Delta x_{i}^{3}}} \right)}} = {{\overset{\hat{}}{\rho}}_{k,i}\Delta x_{i}}} & (26) \end{matrix}$ $\begin{matrix} {{{c_{0}\Delta x_{i + 1}} + {c_{1}{\overset{¯}{x}}_{i + 1}\Delta x_{i + 1}} + {c_{2}\left( {{{\overset{¯}{x}}_{i - 1}^{2}\Delta x_{i - 1}} + {\frac{1}{12}\Delta x_{i - 1}^{3}}} \right)}} = {{\overset{\hat{}}{\rho}}_{k,{i + 1}}\Delta x_{i + 1}}} & (27) \end{matrix}$

Solving eqns. (25) to (27) determines the coefficients c₀ to c₂ which may be used to reconstruct the spatial reconstruction of the mass present in the i'th finite control volume as: {circumflex over (ρ)}_(k,i)*(x)=c₀+c₁x+c₂x², xϵ[x_(i−1/2), x_(i+1/2)]. Combining the polynomial pieces for all the control volumes in the domain give the spatial reconstruction of the mass over the whole domain, {circumflex over (ρ)}_(k,i)*(x), xϵ[x_(1/2), x_(N+1/2)].

Combining the spatially reconstructed mass for the whole domain with the domains of dependence at all internal faces, the mass fluxes between the finite control volumes in Eqn. (4) may be calculated in a new manner. The mass flux over the j'th internal face is calculated by summing the mass over its associated domain of dependence by evaluating the integral

$\begin{matrix} {F_{k,j} = {\frac{A_{j}}{\Delta t_{n}}{{\overset{\hat{}}{\rho}}_{k}^{*}(x)}{dx}}} & (28) \end{matrix}$

where {circumflex over (ρ)}_(k)*(x) is the spatial reconstruction of the mass of the k'th fluid phase over the whole computational domain, composed of all the pieces {circumflex over (ρ)}_(k,i)*(x) for i∈[1, N], which is integrated over the domain of dependence,

_(k,j), and F_(k,j) is the mass flux across the j'th internal face. A_(j) is the cross-sectional area at the position of the j'th internal face.

Higher order polynomials may have regions where they have negative values leading to summation of negative masses. Thus, in an example embodiment, the invention according to the first aspect may further comprise a rescaling of the spatial reconstruction of the mass which may by a limiter function known from [Ref. 2, slide 21] which limits the polynomial to zero or positive values while maintaining the condition defined in eqn. (23):

(x)=θ({circumflex over (ρ)}_(k,i)*(x)−{circumflex over (ρ)}_(k,i) ⁰)+{circumflex over (ρ)}_(k,i) ⁰  (29)

where

$\theta = \frac{{\overset{\hat{}}{\rho}}_{k,i}^{0}}{{\overset{\hat{}}{\rho}}_{k,i}^{0} - m}$

if m<0 or otherwise θ=1, and where

$m = {\min\limits_{x \in {\lbrack{x_{i - {1/2}},x_{i + {1/2}}}\rbrack}}{{{\overset{\hat{}}{\rho}}_{k,i}^{*}(x)}.}}$

Close to the pipe nodes, the domain of dependence may reach outside of the computational domain

. Then, part of the mass flowing through will originate from the outside of the computational domain. Depending on the type of boundary condition that we want to achieve, the following can be done:

In one example embodiment applying an imposed mass flow rate boundary condition, the mass of phase k flowing through internal face j during the remainder of the time step, Δt_(r,k), that was stored in the algorithm to determine the domain of dependence. will be:

F _(in,k) Δt _(r,k)

where F_(in,k) is the imposed mass flow rate. Then, the mass flow rate through the internal face during time step Δt will be determined by integral (28) to which is summed the flow rate during the remainder of the time step, as

$\begin{matrix} {F_{k,j} = {\frac{1}{\Delta t_{n}}\left( {{A_{j}{{\overset{\hat{}}{\rho}}_{k}^{*}(x)}{dx}} + {F_{{in},k}\Delta t_{r,k}}} \right)}} & \left( {28a} \right) \end{matrix}$

where

_(k,j) ∩

denotes the part of the domain of dependence which is within the computational domain.

In another example embodiment applying an imposed pressure boundary condition, the velocity may be applied to decide the mass flow rate, given a user-defined phase split of mass entering the pipe. In turn, the difference between the user-imposed pressure and the pipe's pressure decides the velocity update during the time step, through the pressure gradient term. This may be handled by constantly extrapolating the velocity and applying Equation (17) with ΔC_(k,i)=0, and ū_(k,i) being the velocity in the first or last cell, for the west or east boundary condition, respectively. Equation (17) gives in this case a x_(s,k) ^(n) outside of the computational domain. Then, the mass flow rate through the internal face during time step Δt may be determined by integral (28), in which for the part of the domain of dependence outside the computational domain it is applied a mass defined as {circumflex over (ρ)}(x)=α_(k,in)ρ_(k,in), where the index in denotes the volume fraction imposed at the node and the density corresponding to the imposed pressure.

Numerical schemes which are higher order in space are known to generate spurious oscillations at discontinuities and extrema. On the other hand, numerical schemes which are first order in space do not. Thus, in one embodiment, the rescaling of the spatial reconstruction of the mass further comprise a step to modulate the order of the scheme locally using limiters, based, in each cell, on the value of the solution in the cell and its immediate neighbours. This gives a combined spatially first- and higher-order numerical scheme.

The spatial reconstruction of the mass may be made locally zeroth order by setting θ=0 in eqn. (29), that is, make the reconstruction constant in a cell. A zeroth-order mass reconstruction gives a spatially first-order numerical scheme. Thus, the limiter is defined to prevent the reconstructed polynomial from overshooting or undershooting the neighbouring cell values, by using a limiting condition of the same form as the one limiting negative values, i.e. applying eqn. (29) where, θ is determined as follows:

At the left internal face of the i'th internal control volume, set:

$\theta_{L} = {{❘\frac{{\hat{\rho}}_{k,{i - 1}}^{0} - {\hat{\rho}}_{k,i}^{0}}{{{\hat{\rho}}_{k,i}^{*}\left( x_{i - {1/2}} \right)} - {\hat{\rho}}_{k,i}^{0}}❘}{if}{either}\left\{ \begin{matrix} {{\hat{\rho}}_{k,i}^{0} > {{\hat{\rho}}_{k,{i - 1}}^{0}{and}{{\hat{\rho}}_{k,i}^{*}\left( x_{i - {1/2}} \right)}} < {\hat{\rho}}_{k,{i - 1}}^{0}} \\ {{\hat{\rho}}_{k,i}^{0} < {{\hat{\rho}}_{k,{i - 1}}^{0}{and}{{\hat{\rho}}_{k,i}^{*}\left( x_{i - {1/2}} \right)}} > {\hat{\rho}}_{k,{i - 1}}^{0}} \end{matrix} \right.}$

and at the right internal face, set:

$\theta_{R} = {{❘\frac{{\hat{\rho}}_{k,{i - 1}}^{0} - {\hat{\rho}}_{k,i}^{0}}{{{\hat{\rho}}_{k,i}^{*}\left( x_{i - {1/2}} \right)} - {\hat{\rho}}_{k,i}^{0}}❘}{if}{either}\left\{ \begin{matrix} {{\hat{\rho}}_{k,i}^{0} > {{\hat{\rho}}_{k,{i - 1}}^{0}{and}{{\hat{\rho}}_{k,i}^{*}\left( x_{i - {1/2}} \right)}} < {\hat{\rho}}_{k,{i - 1}}^{0}} \\ {{\hat{\rho}}_{k,i}^{0} < {{\hat{\rho}}_{k,{i - 1}}^{0}{and}{{\hat{\rho}}_{k,i}^{*}\left( x_{i - {1/2}} \right)}} > {\hat{\rho}}_{k,{i - 1}}^{0}} \end{matrix} \right.}$

and then set θ=MIN(θ_(L), θ_(R)), and in addition, at extrema, that is, for cells satisfying the following condition

({circumflex over (ρ)}_(k,i) ⁰>{circumflex over (ρ)}_(k,i−1) ⁰ and {circumflex over (ρ)}_(k,i) ⁰>{circumflex over (ρ)}_(k,i+1) ⁰)or({circumflex over (ρ)}_(k,i) ⁰<{circumflex over (ρ)}_(k,i−1) ⁰, and {circumflex over (ρ)}_(k,i) ⁰<{circumflex over (ρ)}_(k,i+1) ⁰), set θ=0.

Finally, we keep the most restrictive θ (closest to 0) between the presently defined limiter and the limiter for the positivity preserving property in eqn. (29).

In a second aspect, the invention relates to a method for optimising the design of a pipeline-based fluid transportation system for transporting a multiphase fluid flow, wherein the method comprises:

-   -   preparing at least two different designs of the fluid         transportation system,     -   applying the computer implemented method according to the first         aspect of the invention to predict the fluid behaviour in each         of the at least two different designs, and     -   applying the predicted fluid behaviour to determine the         optimised design of the fluid transportation system.

The optimisation of the design of the transportation system may take into consideration one or more factors such as pipeline diameter, pipeline trajectory in the terrain, number of pumps for pressure support, their location and pressure enhancing effect, number of chocking valves, their location and flow volume reducing effect, etc. with the aim to save capital investment and operational costs by identifying the optimum physical dimensions and/or trajectory in the terrain of the transport systems pipes without compromising on fluid behaviour stability and throughput.

In a third aspect, the invention relates to a method for trouble-shooting flow problems during operation of a pipeline-based fluid transportation system for transporting a multiphase fluid flow, wherein the method comprises:

-   -   applying the computer implemented method according to first         aspect of the invention loaded with a computational domain         representative for section of the transport system having flow         problems and with flow characteristic input data of the flow in         the transportation system to predict the effect on the fluid         behaviour in the transport system from possible mitigation         actions, and     -   applying the predicted fluid behaviours to determine which         mitigation action which is to be physically implemented on the         transport system having flow problems.

The mitigation actions may be regulating he flow volumes in the transport system, topside choking, gas lift, and others.

In a fourth aspect, the invention relates to a computer program, comprising processing instructions which causes a computer to perform the method according to any of the first to the third aspects of the invention when the instructions are executed by a processing device in the computer.

In a fifth aspect, the invention relates to a computer, comprising a processing device and a computer memory, the computer memory is storing a computer program as set forth in the fourth aspect.

LIST OF FIGURES

FIG. 1 illustrates schematically an example embodiment of a pipeline system for transporting processed fluids in oil and gas extraction.

FIGS. 2 a ) and 2 b) illustrate graphically the limitation of the CFL-condition when estimating mass fluxes out of control volumes.

FIG. 3 a ) is a schematic representation of an example of a spatial reconstruction of the mass of the k'th fluid phase in the finite control volumes using a first order polynomial to represent the mass as linear interpolation lines over each finite volume.

FIG. 3 b ) is a schematic representation of the same spatial reconstruction as shown in FIG. 3 a ) and which also illustrates an example of a domain of dependence stretching in counter-flow direction from the internal face at position x_(j) to a position x_(0,k) ^(n) being inside the i−3'th finite control volume.

FIG. 3 c ) is a schematic representation of the same spatial reconstruction as shown in FIGS. 3 a ) and 3 b) and which also illustrates the mass being present in the domain of dependence,

_(k,j), by the shaded area.

FIG. 4 is a graphical presentation of a comparison of predicted fluid behaviour in a pipeline made with a prior art solver and a solver according to the invention.

FIG. 5 is a graphical presentation of a comparison of simulated liquid level (continuous liquid plus bubbles) in a second comparison test of a stretch from 1000 to 2000 meter of a pipeline.

FIG. 6 is a graphical presentation of a comparison of simulated liquid level (continuous liquid plus bubbles) in a third comparison test of a stretch from 1000 to 2000 meter of a pipeline.

FIG. 7 is a graphical presentation of negative volume fractions caused by the minmod limiter used instead of the PPS, in the third comparison test.

FIG. 8 is a graphical presentation of a comparison of simulated volume fraction of continuous liquid, plotted on the top row, and volume fraction of bubbles, plotted on the bottom row, of a fourth comparison test.

FIG. 9 is a graphical presentation of a predicted flow using a solver not using the numerical scheme according to the invention and which the source term for the hydraulic force is deactivated to give the same velocity for the gas and liquid phase.

FIG. 10 is a graphical presentation of a predicted flow using a solver using the numerical scheme according to the invention and which the source term for the hydraulic force is deactivated to give the same velocity for the gas and liquid phase.

VERIFICATION OF THE INVENTION

The numerical scheme according to the invention enables predicting multiphase flows in pipelines with an improved (lower) numerical diffusivity as compared to prior art numerical models. Simulations of hydrodynamic wave growth leading to plug flow (when the wave crests reach the upper parts of the pipe wall and slows the liquid flow) are particularly sensitive to numerical diffusion. At the same time, plug flows involve relatively rapid moving fluid phases involving relatively high pressure and mass gradients which requires relatively fine mesh sizes and relatively small time steps making such simulations particularly computational heavy.

Comparison Test 1

The fluid behaviour in a 400 m long (linear) pipe having an internal diameter of 0.194 m and inclined 1° downward, and which is supplied with gas corresponding to a superficial velocity of 3.8 m/s and a liquid corresponding to a superficial velocity of 1 m/s at pressure 20 bar is predicted with a prior art solver utilising an implicit numerical scheme in which the mass fluxes are made partially explicit and compared with a similar prediction made with an implicit-explicit solver (implicit in pressure and velocity, explicit in mass and momentum convection) having incorporated the numerical scheme according to the first aspect of the invention. These flow conditions are known to cause wave growth towards the end of the pipe.

By performing the predictions with a relatively coarse grid, the numerical diffusion will dominate the physical models of the CFD-code and thus enabling assessing which meshes on the prior art and the present solver yield similar results by comparing the distance required for onset of the wave build-up as a measure of the accuracy of the numerical schemes.

The prior art solver was run three times with a mesh size (Δx_(i)) equal to 0.5, 1, and 2 timed the pipe diameter D, respectively. The result is shown graphically in FIG. 4 under the header “old solver”. The predictions gave a stratified flow at the downstream end of the pipe and an onset of wave build-up at 200 to 300 metres downstream into the pipe. The solver according to the first aspect of the invention obtained the same onset of the wave build-up when applying a mesh size of 7, 10, and 14 times the pipe diameter. These results are also shown in FIG. 4 under the header “new solver”.

From FIG. 4 we see that a solver according to the invention may use a mesh in the order of 10 times coarser to yield a similar onset of wave growth as a prior art solver, which represents a significant reduction in the computational load. The comparison is rather qualitative, but still gives an order of magnitude of the reduction in numerical diffusion, which gives increased numerical accuracy. Note that this case only involves hydrodynamic wave growth, which is particularly sensitive to numerical diffusion. One cannot necessarily expect such mesh ratios in all types of flows, since for example terrain slugging is not particularly sensitive to numerical diffusion. Still, even in cases dominated by terrain slugging, it might still be important to resolve hydrodynamic slugs, thus setting a minimum requirement on the refinement of the mesh.

Comparison Test 2

This test investigates the numerical load when using an explicit solver having implemented the numerical scheme according to the first aspect of the invention to simulate the fluid behaviour in a 2 km long stretch and duration of 500 seconds of the same pipe fed with the same multiphase flow as described above for comparison test 1. The solver applies a second order polynomial for spatial reconstruction of the mass in the finite control volumes and a second order limiter function to rescale the spatial reconstruction to preserve positivity of the mass over the computational domain.

As a comparison, the same explicit solver is applied without the numerical scheme according to the first aspect of the invention but instead applies a higher-order upwind scheme using the minmod limiter function [Ref. 3, page 110]. A second comparison is also made with the same explicit solver without the numerical scheme according to the first aspect of the invention but a higher-order upwind scheme using the monotonized central-difference (MC) limiter function [Ref. 3, p. 112]. MC limiter is expected to give more accurate results than minmod, at the risk of being numerically unstable. Amongst other, positivity of mass is never assured with higher-order limiters, but minmod is very mild and safe. MC is recommended in the same reference to be “a good default choice for a wide class of problems”.

Without the numerical scheme according to the first aspect of the invention, it is necessary to apply a CFL condition defined with the true numerical velocities (the actual u_(k,j) used in the scheme), set to be CFL<0.8. In addition, the CFL may be violated during or at the end of the time step. This cannot be allowed, and if CFL>1 during or at the end of the time step, the time step is restarted with a smaller Δt. Time step restarts can happen even with the numerical scheme according to the first aspect of the invention, for other reasons. Even though the numerical scheme according to the first aspect of the invention removes the CFL restriction related to mass convection, other physical phenomena which may be implemented in the numerical model may require strict time step restriction to preserve numerical stability (e.g. for surface waves), or less strict restriction to preserve modelling accuracy (e.g. friction of phase change). Thus, a CFL-condition is still applied in these example numerical calculations, but less strictly and only at the beginning of the time step. A local violation of the CFL condition is acceptable—avoiding restriction of the time step for the whole pipe due to one too high velocity—as is also a violation of the CFL condition at the end of the time step, and as such a time step restart can be avoided.

The resulting computational load is shown in table 1. Table 1 compares the time steps applied in the numerical simulations, number of restarts during the simulation, total used CPU time and time used to solve the mass transport and pressure equations (marked as “Time in mass conv.”). The column marked “PPS” is the simulation with the numerical scheme according to the first aspect of the invention. The column marked “NO PPS—minmod” is the first comparison simulation without the numerical scheme according to the first aspect of the invention but applying the minmod limiter function and the column marked “No PPs—MC” is the second comparison simulation without the numerical scheme according to the first aspect of the invention but applying the MC limiter function. The mesh size in all three simulations was 5 times the pipe diameter.

TABLE 1 Comparison of numerical workload PPS No PPS - minmod No PPS - MC Time step, Δt 0.0933 s 0.0971 s 0.0686 s Number of restarts 0 90 1356 Used CPU time 246 s 208 s 328 s Time in mass conv. 48 s 7 s 25 s

FIG. 5 is a graphical presentation of the simulated liquid level (continuous liquid plus bubbles) in the stretch from 1000 to 2000 meter of the pipe for all three simulations. The “left box” marked “PPS” presents the simulated liquid level along the pipe segment obtained from the explicit solver having implemented the numerical scheme according to the first aspect of the invention. The “middle box” marked “No PPS—minmod” presents the simulated liquid level along the pipe segment obtained from the explicit solver using the minmod limiter function without the numerical scheme according to the first aspect of the invention, while the “right box” marked “No PPS—MC” presents the same for the simulation with the explicit solver using the MC limiter function without the numerical scheme according to the first aspect of the invention.

The explicit solver having implemented the numerical scheme according to the first aspect of the invention (PPS) runs with the same average time step as the explicit solver using the minmod limiter function (minmod). This means that the ability to avoid time step restriction due to mass convection is not playing a role in this comparison. Most of the time, however, this ability is useful in case of velocity spikes at the slug/bubble transition, which is an artifact of the underlying physical models. The time step is considerably lower with the explicit solver using the MC limiter function (MC), as well as the number of time step restarts, indicating that there are probably velocity spikes with this scheme.

The numerical scheme according to the first aspect of the invention is relatively computationally heavy, as can be observed in the lower row of Table 1. Despite running with the same time step, the CPU time of the numerical scheme according to the first aspect of the invention (PPS) is 18% higher than with the explicit solver using the minmod limiter function (minmod). The difference in run time is coming from the execution time of the PPS, which is what makes the difference in the row ‘Time in mass conv’ of Table 1. The comparison of the CPU time with MC is however in favor of the PPS, due to the lower average time step, resulting in a higher total number of time steps to solve, as well as to the number of time step restarts, which requires to start over with a new smaller time step.

Comparison Test 3

This test is similar to comparison test 2 above except for applying the same three numerical simulations on a 2000 meters long linear pipe having an internal diameter of 0.290 m and inclined 1° upward, and which is supplied with gas corresponding to a superficial velocity of 1.5 m/s and a liquid corresponding to a superficial velocity of 0.075 m/s at pressure 20 bar. The results are given in table 2 and in FIG. 6 , respectively.

TABLE 2 PPS No PPS, minmod No PPS, MC Time step 0.2160 s 0.1106 s Crashed N restarts 188 273 Total CPU time 94 s 124 s Time in mass 26 s 6 s convection module

In this case, the explicit solver using the MC limiter function (MC) was not able to complete the simulation, due to instability issues. The explicit solver using the minmod limiter function ran with a time step close to twice as small as the explicit solver having implemented the numerical scheme according to the first aspect of the invention. The consequence is that the explicit solver with the PPS in this case run 25% faster, i.e. a total CPU time of 94 s versus 124 s, even though it is spending 20 s more, 26 s versus 6 s, in run time solving the mass transport equations as compared to the No PPS run.

These results indicate that “No PPS—minmod” seems to be more diffusive than PPS, as suggested by the fact that many waves are not reaching the top of the pipe and becoming slugs. In addition, negative masses are observed with minmod as limiter, even though the minmod limiter is a mild higher-order limiter. FIG. 7 shows an example of this in a zoom on the simulated gas volume fraction using the No PPS—minmod scheme (middle plot in FIG. 6 ). With the PPS, no negative mass can be observed at any time step.

Comparison Test 4

This test is a 1000 m pipe of diameter 0.1 m, initialised with single-phase gas at rest. Then, liquid is injected at the left node at a superficial velocity of 1 m/s. Both liquid and gas are incompressible. We expect to see a liquid front between single-phase gas and single-phase liquid propagate at a velocity of 1 m/s. The same numerical schemes as in tests 2 and 3 are applied (PPS, No PPS—minmod and No PPS—MC, referring to a solver having implemented the numerical scheme according to the first aspect of the invention, a solver using the minmod limiter function without the numerical scheme according to the first aspect of the invention, and a solver using the MC limiter function without the numerical scheme according to the first aspect of the invention, respectively).

After 900 s, the results are plotted in FIG. 8 . The volume fraction of continuous liquid is plotted on the top row, and the volume fraction of bubbles is plotted on the bottom row. With the PPS applied, there are no volume fractions below 0 or above 1, while both minmod and MC limiters cause negative volume fractions (the sum of all volume fractions−continuous gas and liquid, bubbles and droplets at a given position in the pipe are always equal to 1, plus or minus the convergence criterion-equal to 1e-3 in the present case). This is a well-known behaviour of higher-order limiters to cause oscillations at discontinuities if they do not degenerate to first order fast enough. Minmod being more cautious than MC, it causes milder oscillations. The result of oscillations is that negative volume fractions are predicted. We can see here that the PPS keeps all volume fractions positive.

Comparison Test 5

For this test, the source term for the hydraulic force (pressure gradient due to a gradient in liquid level) is deactivated. Thus, if the gas and liquid velocities are exactly equal, the liquid-gas interface is expected to be advected as is at the same velocity. It is run with a solver having implemented the numerical scheme according to the first aspect of the invention, with a mass reconstruction using second order polynomials. A 300 m flat pipe is meshed with 318 cells, with an initial velocity of 5.57 m/s both in gas and in liquid, and an initial crenel-shaped gas-liquid interface as shown in FIG. 9 , left plot. The liquid and gas densities are respectively 818.7 kg/m³ and 64 kg/m³. The inlet boundary condition is defined to inject the phases with exactly the right mass rate to keep the same volume fractions and velocities. In FIG. 9 , right plot, the result after advection during 30 s is plotted. The crenel-shaped interface is advected and somewhat similar to the initial shape, but spurious oscillations at both discontinuities have appeared. In FIG. 10 , the result is plotted after running the test with a solver having implemented the numerical scheme according to the first aspect of the invention, with a mass reconstruction using second order polynomials, and a limiter to avoid spurious oscillations. On the right plot, the spurious oscillations have disappeared. Numerical diffusion has smoothed the discontinuities to some extent, as expected. The crenel shape is otherwise advected as is. The small waves at ca. 150 m are small disturbances caused by the inlet node at the start of the case, which have been advected.

REFERENCES

-   1 Veersteg, H. K. and Malalasekera, W., sec. Ed. (2007), “An     Introduction to Computational Fluid Dynamics”, Pearson Education     Limited, ISBN: 978-0-13-127498-3 -   2 Chi-Wang Shu, “Positivity-preserving high order schemes for     convection dominated equations”, presentation held at Brown     University and retrievable on the internet:     https://cfd.ku.edu/JRV/Shu.pdf -   3 Randall J. Leveque, Finite Volume Methods for Hyperbolic Problems,     Cambridge University Press, 2002, ISBN 978-O-521-00924-9 -   4 M. Griebel and M. Klitz, “CLSVOF as a fast and mass-conserving     extension of the level-set method for the simulation of two-phase     flow problems”, NUMERICAL HEAT TRANSFER, PART B 2017, VOL. 71, NO.     1, 1-36 http://dx.doi.org/10.1080/10407790.2016.1244400 

1. A computer implemented method for predicting fluid behaviour of a multiphase flow in a pipeline-based transport system where the flow contains a number of k, where k is a positive integer, fluid phases, wherein the method comprises: applying a one-dimensional (1D) computational fluid dynamic (CFD) model describing the geometry of a section of interest of the pipeline-based transport system and the multiphase flow flowing therein, solving the 1D CFD model to simulate the fluid behaviour of the multiphase flow in the section of interest of the pipeline-based transport system, and describing the determined fluid behaviour by presenting the simulation results as macroscopic fluid properties such as flow velocity, pressure, density, temperature, etc., and the 1D CFD model applies a finite volume method to solve the model, wherein the geometry of the section of interest of the pipeline-based transport system is defined as a computational domain extending along an axis represented by the cartesian coordinate x and being divided into a set of N, where N is a positive integer, non-overlapping finite control volumes separated by an internal face between adjacent finite control volumes, characterised in that the one-dimensional computational fluid dynamic model is adapted to estimate the mass flow of a k'th fluid phase out of a i'th finite control volume during a n'th time step by applying a polynomial to spatially reconstruct the mass, {circumflex over (ρ)}_(k)(x), of the k'th fluid phase being present in each of the N finite control volumes of the computational domain, and then for each j'th internal face, where jϵ1/2, . . . , N+1/2, of the computational domain, execute the following steps i) and ii): i) reconstruct the flow velocity, u_(k)(x), of the k'th fluid phase as a function of position x and apply the reconstruction to determine a domain of dependence,

_(k,j), representing the distance the k'th fluid phase has travelled during the n'th time step, and ii) sum the spatially reconstructed mass being present in the domain of dependence,

_(k,j), and assume the summarised mass passes through the j'th internal face over the n'th time step, into the i'th finite control volume when u_(k)(x_(j))<0 or into the i+1'th finite control volume when u_(k)(x_(j))>0, where u_(k)(x_(j)) is the flow velocity at the j'th internal face.
 2. The computer implemented method according to claim 1, wherein the CFD-model applies an explicit numerical solution scheme.
 3. The computer implemented method according to claim 1, wherein the method further comprises determining the domain of dependence for the i'th finite control volume,

_(k,j) by executing the following steps i) to vii): i) if u_(k,j)>0: set an upwind cell index i=j−1/2 and a direction index s=−1, if u_(k,j)<0: set an upwind cell index i=j+1/2 and a direction index s=1, ii) set a cell counter m=0, iii) set Δt_(r) equal to the full time step Δt, iv) if |r_(k,i+sm)−1|<2√{square root over (3ϵ)}: ${{{set}\Delta t_{c,k,{i + {sm}}}} = \frac{\Delta x_{i + {sm}}}{{\overset{¯}{u}}_{k,{i + {sm}}}}},$ otherwise: ${{{set}\Delta t_{c,k,{i + {sm}}}} = {\frac{\Delta x_{i + {sm}}}{{\overset{¯}{u}}_{k,{i + {sm}}}}\left( {\frac{r_{k,{i + {sm}}} + 1}{2}\frac{\ln\left( r_{k,{i + {sm}}} \right)}{\left( {r_{k,{i + {sm}}} - 1} \right)}} \right)}},$ where ${r_{k,{i + {sm}}} = \frac{u_{k,{i + {sm} + {1/2}}}}{u_{k,{i + {sm} - {1/2}}}}},$  and ϵ is a positive real number obtained from executing a Fortran EPSILON computer implemented function, v) if Δt_(c,k,i+sm)>Δt_(r): go to step vi), or else if i+s(m+1)<1 or i+s(m+1)>N: set Δt_(r)=Δt_(r)−Δt_(c,k,i+sm), and terminate the procedure, or else: set m=m+1 and Δt_(r)=Δt_(r)−Δt_(c,k,i+sm), and go to step iv), vi) if |ΔC_(k,i+sm)|<√{square root over (6ϵ)}: determine a characteristic starting position, x_(0,k) ^(n), by solving the following eqn.; ${x_{0,k}^{n} = {{\overset{\_}{x}}_{i + {sm}} + {e^{{- \Delta}C_{k,{i + {sm}}}}\left( {x_{j + {sm}} - {\overset{\_}{x}}_{i + {sm}}} \right)} - {{\overset{\_}{u}}_{k,{i + {sm}}}\Delta{t_{r}\left( {1 - {\frac{1}{2}\Delta C_{k,{i + {sm}}}}} \right)}}}},$  and go to step vii), Otherwise: determine a characteristic starting position, x_(0,k) ^(n), by solving the following eqn.; ${x_{0,k}^{n} = {{\overset{\_}{x}}_{i + {sm}} + {e^{{- \Delta}C_{k,{i + {sm}}}}\left( {x_{j + {sm}} - {\overset{\_}{x}}_{i + {sm}}} \right)} - {{\overset{\_}{u}}_{k,{i + {sm}}}\Delta{t_{r}\left( \frac{\left( {1 - e^{{- \Delta}C_{k,{i + {sm}}}}} \right)}{\Delta C_{k,{i + {sm}}}} \right)}}}},$  and go to step vii), where ${{\overset{\_}{x}}_{i + {sm}} = \frac{x_{i + {sm} + {1/2}} + x_{i + {sm} - {1/2}}}{2}},{{\overset{\_}{u}}_{k,{i + {sm}}} = \frac{u_{k,{i + {sm} - {1/2}}} + u_{k,{i + {sm} + {1/2}}}}{2}},$ Δx_(i + sm) = x_(i + sm + 1/2) − x_(i + sm − 1/2), Δu_(k, i + sm) = u_(k, i + sm + 1/2) − u_(k, i + sm − 1/2), ${{\Delta C_{k,{i + {sm}}}} = {\Delta u_{k,{i + {sm}}}\frac{\Delta t}{\Delta x_{i + {sm}}}}},$  and vii) if u_(k,j)>0: set

_(k,j)=[x_(0,k) ^(n), x_(j)], or if u_(k,j)<0: set

_(k,j)=[x_(j), x_(0,k) ^(n)].
 4. The computer implemented method according to claim 3, wherein the spatial reconstruction over the whole domain of the mass, {circumflex over (ρ)}_(k)*(x), of the k'th fluid phase being present in each of the N finite control volumes of the computational domain is obtained by: applying in each cell a polynomial of even degree: ${{{\hat{\rho}}_{k,i}^{*}(x)} = {\sum\limits_{\alpha = 0}^{\beta}{c_{\alpha}x^{\alpha}}}},{\beta \in \left\lbrack {0,2,4,\ldots} \right\rbrack}$ and further by: for each i'th finite control volume, i=1 to N: define a set of coefficients, c_(α), through the condition: ${{\overset{x_{i + p + {1/2}}}{\int\limits_{x_{i + p - {1/2}}}}{{{\hat{\rho}}_{k,i}^{*}(x)}{dx}}} = {{\hat{\rho}}_{k,{i + p}}\Delta x_{i + p}}},$ where p∈[−β/2, β/2] is an integer number, and solve the resulting system of equations for the coefficients c₀, c₁, . . . , c_(β) to reconstruct the spatial reconstruction of the mass, {circumflex over (ρ)}_(k,i)*(x), present in the i'th finite control volume as: {circumflex over (ρ)}_(k,i)*(x)=c₀+c₁x+c₂x²+ . . . +c_(β)x^(β), xϵ[x_(i−1/2), x_(i+1/2].)
 5. The computer implemented method according to claim 4, wherein the even degree polynomial is a second order polynomial, c₀+c₁x+c₂x², and using that: ${{\overset{x_{i + {1/2}}}{\int\limits_{x_{i - {1/2}}}}c_{0}} + {c_{1}x} + {c_{2}x^{2}{dx}}} = {{c_{0}\Delta x_{i}} + {c_{1}{\overset{\_}{x}}_{i}\Delta x_{i}} + {c_{2}\left( {{{\overset{\_}{x}}_{i}^{2}\Delta x_{i}} + {\frac{1}{12}\Delta x_{i}^{3}}} \right)}}$ to define relations for the i−1'th, the i'th, and the i+1'th finite control volumes, respectively: ${{c_{0}\Delta x_{i - 1}} + {c_{1}{\overset{\_}{x}}_{i - 1}\Delta x_{i - 1}} + {c_{2}\left( {{{\overset{\_}{x}}_{i - 1}^{2}\Delta x_{i - 1}} + {\frac{1}{12}\Delta x_{i - 1}^{3}}} \right)}} = {\hat{\rho}}_{k,{i - 1}}$ ${{c_{0}\Delta x_{i}} + {c_{1}{\overset{\_}{x}}_{i}\Delta x_{i}} + {c_{2}\left( {{{\overset{\_}{x}}_{i}^{2}\Delta x_{i}} + {\frac{1}{12}\Delta x_{i}^{3}}} \right)}} = {\hat{\rho}}_{k,i}$ ${{c_{0}\Delta x_{i + 1}} + {c_{1}{\overset{\_}{x}}_{i + 1}\Delta x_{i + 1}} + {c_{2}\left( {{{\overset{\_}{x}}_{i + 1}^{2}\Delta x_{i + 1}} + {\frac{1}{12}\Delta x_{i + 1}^{3}}} \right)}} = {\hat{\rho}}_{k,{i + 1}}$ and solving the three second order polynomials to determine the coefficients c₀, c₁ and c₂, and then reconstruct the spatial reconstruction of the mass, {circumflex over (ρ)}_(k,i)*(x), present in the i'th finite control volume as: {circumflex over (ρ)}_(k,i)*(x)=c₀+c₁x+c₂x², xϵ[x_(i−1/2), x_(i+1/2)].
 6. The computer implemented method according to claim 4, wherein the spatially reconstructed mass is applied to estimate a mass flux, F_(k,j), across the j'th internal face by: $F_{k,j} = {\frac{A_{j}}{\Delta t_{n}}{\int_{{\mathbb{D}}_{k,j}}{{{\hat{\rho}}_{k}^{*}(x)}{dx}}}}$ where {circumflex over (ρ)}_(k)*(x) is the reconstruction of the mass of the k'th fluid phase over the whole computational domain, composed of the polynomials {circumflex over (ρ)}_(k,i)*(x) from all cells, integrated over the j'th domain of dependence,

_(k,j), Δt_(n) is the n'th time step, and A_(j) is the cross-sectional area at the position of the j'th internal face.
 7. The computer implemented method according to claim 6, wherein the method further comprises, when applying an imposed mass flow rate, F_(in,k), into the computational domain, that for each internal face j having a domain of dependence,

_(k,j) with a starting point, x_(0,k) ^(n), being outside of the computational domain

, the mass flow rate through internal face j is determined as: $F_{k,j} = {\frac{1}{\Delta t_{n}}\left( {{A_{j}{\int\limits_{{\mathbb{D}}_{k,j}\bigcap{\mathbb{P}}}{{{\hat{p}}_{k}^{*}(x)}{dx}}}} + {F_{{in},k}\Delta t_{r,k}}} \right)}$ wherein

_(k,j)∩

denotes the part of the domain of dependence which is within the computational domain.
 8. The computer implemented method according to claim 6, wherein the method further comprises, when applying an imposed pressure boundary condition, that for each internal face j having a domain of dependence,

_(k,j) with a starting point, x_(0,k) ^(n), being outside of the computational domain

, the mass flow rate through internal face j is determined by extrapolating the velocity and applying: $x_{s,k}^{n} = {{\overset{\_}{x}}_{i} + {e^{{- \Delta}C_{k,i}}\left( {x_{j} - {\overset{\_}{x}}_{i}} \right)} - {{\overset{\_}{u}}_{k,i}\Delta t\left\{ \begin{matrix} {1 - {\frac{1}{2}\Delta C_{k,i}}} & {{{if}{❘{\Delta C_{k \cdot i}}❘}} < \sqrt{6\epsilon}} \\ \frac{\left( {1 - e^{{- \Delta}C_{k,i}}} \right)}{\Delta C_{k,i}} & {otherwise} \end{matrix} \right.}}$ with ΔC_(k,i)=0, and ū_(k,i) being the velocity in the first or last finite control volume of the computational domain at the west or east boundary condition, respectively to determine an updated starting position, x_(s,k) ^(n), being outside of the computational domain, and then determining the mass flow rate through the internal face j during time step Δt_(n) by: $F_{k,j} = {\frac{A_{j}}{\Delta t_{n}}{\int\limits_{{\mathbb{D}}_{k,j}}{{{\hat{p}}_{k}^{*}(x)}{dx}}}}$ in which for the part of the domain of dependence outside the computational domain, it is applied a mass defined as {circumflex over (ρ)}_(k)*(x)=α_(k,in)ρ_(k,in), where α_(k,in) is the volume fraction of phase k imposed at the internal face j and ρ_(k,in) is a density corresponding to the imposed pressure.
 9. The computer implemented method according to claim 4, wherein the method further comprises, for each i'th finite control volume, i=1 to N, a rescaling of the polynomial {circumflex over (ρ)}_(k,in)(x) to preserve positivity by:

(x)=θ({circumflex over (ρ)}_(k,i)*(x)−{circumflex over (ρ)}_(k,i) ⁰)+{circumflex over (ρ)}_(k,i) ⁰ where

is the rescaled polynomial for the i'th finite control volume, {circumflex over (ρ)}_(k,i) ⁰ is the mass present in the i'th finite control volume at the beginning of the n'th time step, and ${\theta = {{{\frac{{\hat{\rho}}_{k,i}^{0}}{{\hat{\rho}}_{k,i}^{0} - m}{if}m} < {0{or}\theta}} = {{1{if}m} \geq 0}}},{where}$ ${m = {\min\limits_{x \in {\lbrack{x_{i - {1/2}},x_{i + {1/2}}}\rbrack}}{{\hat{\rho}}_{k,i}^{*}(x)}}},$ and then applying the rescaled polynomials

for the i'th finite control volumes, i=1 to N, in the spatial reconstruction of the mass.
 10. The computer implemented method according to claim 4, wherein the method further comprises, for each i'th finite control volume, i=1 to N, a rescaling of the polynomial {circumflex over (ρ)}_(k,i)*(x) to avoid spurious oscillations at discontinuities and extrema by:

(x)=θ({circumflex over (ρ)}_(k,i)*(x)−{circumflex over (ρ)}_(k,i) ⁰)+{circumflex over (ρ)}_(k,i) ⁰ where

is the rescaled polynomial for the i'th finite control volume, {circumflex over (ρ)}_(k,i) ⁰ is the mass present in the i'th finite control volume at the beginning of the n'th time step, and where θ = 0, if(ρ̂_(k, i)⁰ > ρ̂_(k, i − 1)⁰andρ̂_(k, i)⁰ > ρ̂_(k, i + 1)⁰)or (ρ̂_(k, i)⁰ < ρ̂_(k, i − 1)⁰andρ̂_(k, i)⁰ < ρ̂_(k, i + 1)⁰), or θ = MIN(θ_(L), θ_(R)), where $\theta_{L} = {{❘\frac{{\hat{\rho}}_{k,{i - 1}}^{0} - {\hat{\rho}}_{k,i}^{0}}{{{\hat{\rho}}_{k,i}^{*}\left( x_{i - {1/2}} \right)} - {\hat{\rho}}_{k,i}^{0}}❘}{if}{either}}$ $\left\{ {\begin{matrix} {{\hat{\rho}}_{k,i}^{0} > {{\hat{\rho}}_{k,{i - 1}}^{0}{and}{{\hat{\rho}}_{k,i}^{*}\left( x_{i - {1/2}} \right)}} < {\hat{\rho}}_{k,{i - 1}}^{0}} \\ {{\hat{\rho}}_{k,i}^{0} < {{\hat{\rho}}_{k,{i - 1}}^{0}{and}{{\hat{\rho}}_{k,i}^{*}\left( x_{i - {1/2}} \right)}} > {\hat{\rho}}_{k,{i - 1}}^{0}} \end{matrix}.{and}} \right.$ $\theta_{R} = {{❘\frac{{\hat{\rho}}_{k,{i + 1}}^{0} - {\hat{\rho}}_{k,i}^{0}}{{{\hat{\rho}}_{k,i}^{*}\left( x_{i + {1/2}} \right)} - {\hat{\rho}}_{k,i}^{0}}❘}{if}{either}}$ $\left\{ {\begin{matrix} {{\hat{\rho}}_{k,i}^{0} > {{\hat{\rho}}_{k,{i + 1}}^{0}{and}{\hat{\rho}}_{k,i}^{*}\left( x_{i + {1/2}} \right)} < {\hat{\rho}}_{k,{i + 1}}^{0}} \\ {{\hat{\rho}}_{k,i}^{0} < {{\hat{\rho}}_{k,{i + 1}}^{0}{and}{\hat{\rho}}_{k,i}^{*}\left( x_{i + {1/2}} \right)} > {\hat{\rho}}_{k,{i + 1}}^{0}} \end{matrix},} \right.$ and then applying the rescaled polynomials

for the i'th finite control volumes, i=1 to N, in the spatial reconstruction of the mass.
 11. The computer implemented method according to claim 9, wherein the spatially reconstructed mass is applied to estimate a mass flux, F_(k,j), across the j'th internal face by: $F_{k,j} = {\frac{A_{j}}{\Delta t_{n}}{\int_{{\mathbb{D}}_{k,j}}{(x){dx}}}}$ where

(x) is the reconstruction of the mass of the k'th fluid phase over the whole computational domain, composed of the polynomials

(x) from all cells, which may have been rescaled, integrated over the j'th domain of dependence,

_(k,j), Δt_(n) is the n'th time step, and A_(j) is the cross-sectional area at the position of the j'th internal face.
 12. A method for optimising the design of a pipeline-based fluid transportation system for transporting a multiphase fluid flow, wherein the method comprises: preparing at least two different designs of the fluid transportation system, applying the computer implemented method according to claim 1 to predict the fluid behaviour in each of the at least two different designs, and applying the predicted fluid behaviour to determine the optimised design of the fluid transportation system.
 13. The method according to claim 12, wherein the optimisation of the design of the transportation system assesses the effect on the fluid behaviour of varying one or more factors chosen from; pipeline diameter, pipeline trajectory in the terrain, number of pumps for pressure support, their location and pressure enhancing effect, and number of chocking valves, their location and flow volume reducing effect with the aim to save capital investment and operational costs by identifying the optimum physical dimensions and/or trajectory in the terrain of the transport systems pipes without compromising on fluid behaviour stability and throughput.
 14. A method for trouble-shooting flow problems during operation of a pipeline-based fluid transportation system for transporting a multiphase fluid flow, wherein the method comprises: applying the computer implemented method according to claim 1 loaded with a computational domain representative for section of the transport system having flow problems and with flow characteristic input data of the flow in the transportation system to predict the effect on the fluid behaviour in the transport system from possible mitigation actions, and applying the predicted fluid behaviours to determine which mitigation action which is to be physically implemented on the transport system having flow problems.
 15. A computer program, comprising processing instructions which causes a computer to perform the method according to claim 1 when the instructions are executed by a processing device in the computer.
 16. A computer, comprising a processing device and a computer memory, the computer memory is storing a computer program as set forth in claim
 15. 